home *** CD-ROM | disk | FTP | other *** search
/ BBS in a Box 3 / BBS in a box - Trilogy III.iso / Files / Prog / D-G / GemsI / Src / Sturm / util.c < prev   
Encoding:
C/C++ Source or Header  |  1992-06-16  |  1.7 KB  |  122 lines  |  [TEXT/MPS ]
  1.  
  2. /*
  3.  * util.c
  4.  *
  5.  *    some utlity functions for root polishing and evaluating
  6.  * polynomials.
  7.  */
  8. #include <math.h>
  9. #include <stdio.h>
  10. #include "solve.h"
  11.  
  12. /*
  13.  * evalpoly
  14.  *
  15.  *    evaluate polynomial defined in coef returning its value.
  16.  */
  17. double
  18. evalpoly (ord, coef, x)
  19.     int        ord;
  20.     double    *coef, x;
  21. {
  22.     double    *fp, f;
  23.  
  24.     fp = &coef[ord];
  25.     f = *fp;
  26.  
  27.     for (fp--; fp >= coef; fp--)
  28.     f = x * f + *fp;
  29.  
  30.     return(f);
  31. }
  32.  
  33.  
  34. /*
  35.  * modrf
  36.  *
  37.  *    uses the modified regula-falsi method to evaluate the root
  38.  * in interval [a,b] of the polynomial described in coef. The
  39.  * root is returned is returned in *val. The routine returns zero
  40.  * if it can't converge.
  41.  */
  42. int
  43. modrf(ord, coef, a, b, val)
  44.     int        ord;
  45.     double    *coef;
  46.     double    a, b, *val;
  47. {
  48.     int        its;
  49.     double    fa, fb, x, lx, fx, lfx;
  50.     double    *fp, *scoef, *ecoef;
  51.  
  52.     scoef = coef;
  53.     ecoef = &coef[ord];
  54.  
  55.     fb = fa = *ecoef;
  56.     for (fp = ecoef - 1; fp >= scoef; fp--) {
  57.         fa = a * fa + *fp;
  58.         fb = b * fb + *fp;
  59.     }
  60.  
  61.     /*
  62.      * if there is no sign difference the method won't work
  63.      */
  64.     if (fa * fb > 0.0)
  65.         return(0);
  66.  
  67.     if (fabs(fa) < RELERROR) {
  68.         *val = a;
  69.         return(1);
  70.     }
  71.  
  72.     if (fabs(fb) < RELERROR) {
  73.         *val = b;
  74.         return(1);
  75.     }
  76.  
  77.     lfx = fa;
  78.     lx = a;
  79.  
  80.  
  81.     for (its = 0; its < MAXIT; its++) {
  82.  
  83.         x = (fb * a - fa * b) / (fb - fa);
  84.  
  85.         fx = *ecoef;
  86.         for (fp = ecoef - 1; fp >= scoef; fp--)
  87.                 fx = x * fx + *fp;
  88.  
  89.         if (fabs(x) > RELERROR) {
  90.                 if (fabs(fx / x) < RELERROR) {
  91.                     *val = x;
  92.                     return(1);
  93.                 }
  94.         } else if (fabs(fx) < RELERROR) {
  95.                 *val = x;
  96.                 return(1);
  97.         }
  98.  
  99.         if ((fa * fx) < 0) {
  100.                 b = x;
  101.                 fb = fx;
  102.                 if ((lfx * fx) > 0)
  103.                     fa /= 2;
  104.         } else {
  105.                 a = x;
  106.                 fa = fx;
  107.                 if ((lfx * fx) > 0)
  108.                     fb /= 2;
  109.         }
  110.  
  111.         lx = x;
  112.         lfx = fx;
  113.     }
  114.  
  115.     fprintf(stderr, "modrf overflow %f %f %f\n", a, b, fx);
  116.  
  117.     return(0);
  118. }
  119.     
  120.  
  121.  
  122.